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0.1  INTRODUCTION 


This  is  a  final  report  for  grant  (F49620-93- 1-0424)  titled  “Characterization  of  Seismic 
Sources  Using  Empirical  Green’s  Functions”  under  the  Parent  Contract  titled  “Seismic 
Wave  Radiation,  Propagation  and  Event  Location  in  Laterally  Heterogeneous  Media”.  The 
goal  of  this  work  was  to  investigate  the  use  of  the  empirical  Green’s  function  method  for 
seismic  source  characterization  and  discrimination  in  nuclear  monitoring.  The  motivation 
underlying  this  investigation  is  to  support  the  development  of  techniques  capable  of  decom¬ 
posing  complex  seismic  events,  detecting  anomalous  events,  and  discriminating  underground 
nuclear  explosions  from  other  seismic  events. 

Seismic  monitoring  of  the  Comprehensive  Test  Ban  Treaty  (CTBT)  entails  using  teleseismi- 
cally  and  regionally  recorded  ground  motions  generated  by  a  seismic  source  to  estimate  the 
location,  focal  mechanism,  and  time  history  of  energy  release  of  the  source.  These  source 
parameters  allow  the  discrimination  between  underground  nuclear  explosions  and  natural  or 
other  man-made  events.  The  fundamental  difficulty  implicit  in  this  task  is  that  seismic  data 
depend  as  much  on  the  properties  of  the  earth  through  which  the  energy  has  propagated  as 
they  do  on  the  source  parameters.  Thus,  seismologists  attempt  either  to  isolate  features  of 
seismic  data  that  are  relatively  less  sensitive  to  the  propagation  path,  or  to  calibrate  path 
effects  empirically  and  with  numerical  models. 

A  potentially  useful  way  to  deal  with  this  problem  is  to  use  the  recordings  from  multiple 
sources  to  infer  differences  between  their  source  parameters.  That  is,  if  two  or  more  seismic 
events  have  sufficiently  similar  locations,  the  propagation  effects  to  a  given  receiver  will  be 
similar  and  differences  between  the  recorded  seismograms  will  reflect,  primarily,  differences 
between  the  source  parameters. 

One  method  that  adopts  this  approach  is  the  empirical  Green’s  function  method  for  estimat¬ 
ing  the  source  time  function  of  one  event  relative  to  another  (Hartzell,  1978;  Mueller,  1985). 
A  second  method  is  relative  event  location  using  differential  arrival  times  between  event 
pairs,  as  measured  by  the  cross-correlation  of  waveforms  (Poupinet  et  al.,  1984;  Phillips 
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et  al .,  1992).  M.I.T.  has  tested  the  usefulness  of  empirical  Green’s  function  and  relative 
event  location  methods  for  locating  and  identifying  events  in  the  context  of  the  CTBT  in 
previous  work  supported  by  the  Air  Force  (Toksoz  et  al. ,  1993;  Riviere-Barbier  et  al.,  1993; 
Rodi  et  al.,  1994;  Li  et  al.,  1994,  1995, 1996). 

The  EGF  method  attributes  waveform  differences  to  variations  in  the  source  time  functions 
of  two  events.  Thus,  as  it  is  usually  applied,  it  is  required  that  the  two  sources  be  located 
close  together  so  that  none  of  the  waveform  differences  are  due  to  location  or  path  effects. 

In  the  relative  event  location  method,  it  is  assumed  that  seismograms  from  two  closely 
spaced  events  will  have  a  high  degree  of  correlation,  and  that  the  lag  of  their  cross  correlation 
function,  at  a  given  receiver,  reflects  the  difference  between  the  traveltimes  to  the  receiver. 
In  fact,  this  will  be  the  case  only  when  the  source  time  functions  of  the  events  are  the  same 
and  when  scattering  near  the  source  has  the  same  effects  on  their  seismograms.  Waveform 
similarity  between  closely  spaced  events  of  similar  mechanism  is  routinely  observed,  but  the 
degree  of  similarity  appears  to  degrade  with  event  separation.  Israelsson  (1990)  reported 
this  to  be  the  case  for  a  cluster  of  industrial  explosions  he  studied,  while  Riviere-Barbier  and 
Grant  (1991)  employed  this  property  of  waveform  correlation  in  a  cluster-analysis  approach 
to  event  discrimination  .  Results  of  Toksoz  et  al.  (1991a) ,  who  analyzed  waveform  similarity 
as  a  function  of  frequency  and  as  a  function  of  receiver  separation  (for  a  common  earthquake 
source),  suggest  that  scattering  differences  is  the  main  reason  why  waveform  similarity 
degrades  with  event  separation.  Therefore,  one  expects  relative  event  location  to  become 
more  difficult  as  the  distance  between  events  increases. 

In  an  effort  to  develop  seismological  techniques  that  can  be  used  to  discern  natural  earth¬ 
quakes  from  induced  events  at  regional  distances  we  can  consider  the  application  of  either 
EGFs  or  relative  event  location  methods.  However,  as  we  discussed  above,  the  degree  of 
correlation  can  be  shown  to  decrease  with  increasing  distance  between  events  and  is  related 
to  the  degree  of  scattering  along  the  path.  The  problem  of  obtaining  correlations  between 
waveforms  that  represent  true  variations  in  either  location  or  source  time  function  (rather 
than  scattering  effects)  presents  a  significant  challenge  to  the  practical  application  of  either 
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technique. 


In  the  next  section  of  this  report,  we  present  the  results  of  the  application  of  relative  event 
location  to  a  cluster  of  natural  events  as  recorded  by  a  regional  seismic  network  in  the 
Larderello  Geothermal  region  in  Western  Italy.  We  will  show  the  entire  analysis  of  this 
cluster  and  interpret  the  results  with  respect  to  the  known  geology.  Upon  conclusion,  we 
will  review  the  results  of  the  waveform  correlation  for  this  case  study  and  identify  where 
the  practical  limitations  arose  in  using  this  particular  correlation  technique  in  a  regional 
setting. 
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0.2  CASE  STUDY:  LARDERELLO,  ITALY  1993 
MICROEARTHQUAKE  SWARM 


On  March  20, 1993,  seismometers  connected  to  the  Larderello  seismic  network  began  record¬ 
ing  a  series  of  microearthquakes  near  the  Larderello  geothermal  field  in  western  Tuscany, 
Italy  (Figure  0-1.  ENEL  (the  Italian  Electric  Power  Company)  produced  the  initial  esti¬ 
mates  of  the  hypocentral  locations,  the  individual  event  origin  times,  and  estimates  of  their 
magnitudes  (Figures  0-2  and  0-3)  (ENEL-Unita  Nazionale  Geotermica,  1993).  The  Earth 
Resources  Laboratory  (ERL)  at  the  Massachusetts  Institute  of  Technology  (MIT)  extended 
the  analysis  of  this  cluster  to  accomplish  three  objectives;  1)  to  more  precisely  determine  the 
hypocentral  locations  using  multiple  relative  event  location  techniques,  such  that  the  cluster 
geometry  and  overall  location  might  suggest  the  geologic  conditions  and  structures  causing 
the  cluster,  2)  to  determine  focal  mechanisms  for  a  portion  of  the  earthquakes  to  ascertain 
the  local  state  of  stress  and  supplement  the  information  from  the  relocation  inversion,  and 
3)  to  interpret  the  results  from  1  and  2  in  light  of  available  geologic  and  geophysical  data. 

0.2.1  The  Method  of  Multiple  Relative  Event  Location  by  Differential 
Travel  Times 

Applicability  of  the  Method 

The  technique  of  locating  earthquake  hypocenters  using  multiple  relative  event  techniques 
improves  upon  traditional  methods  of  single  event  hypocenter  location  by  placing  additional 
constraints  on  the  results  of  the  inversion.  The  constraints,  which  are  added  to  the  inversion 
in  the  form  of  differential  travel  times,  supplement  the  usual  P  and  S  arrival  times.  The 
incorporation  of  these  additional  constraints  is  suggested  by  the  similarity  of  waveforms 
between  related  events.  The  fundamental  premise  is  that  events  with  identical  waveforms 
must  be  collocated  (  or  very  nearly  so)  (Geller  and  Mueller,  1980;  Poupinet  et  al.,  1984;  Ito, 
1985;  Thornjarnardottir  and  Pechmann,  1987;  Ito,  1990;  Israelsson,  1990;  Pechmann  and 
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Thorbjarnardottir,  1990;  Deichmann  and  Garcia-Fernandez,  1992).  We  can  use  this  simi¬ 
larity  to  improve  the  estimates  of  event  source  parameters  (i.e.  hypocenter  coordinates  and 
origin  times)  if  we  can  identify  and  compensate  for  events  whose  waveforms  are  dissimilar. 
To  do  this  we  must  first  consider  the  cause  of  waveform  variations. 

Fremont  and  Malone  (1987)  state  that  variations  in  waveforms  can  be  attributed  to  four 
main  causes;  fluctuations  in  the  source  time  functions,  changes  in  focal  mechanisms,  changes 
in  the  properties  of  the  media,  and  changes  in  hypocentral  location.  In  order  to  use  waveform 
similarity  to  improve  the  locations,  the  method  of  differential  travel  times  requires  that  the 
set  of  seismic  events  to  be  located  are  related  to  each  other  in  both  spatial  and  temporal 
aspects.  Placing  this  constraint  on  the  events  controls,  to  a  degree,  one  of  the  possible 
sources  of  waveform  variation,  specifically  variations  introduced  by  changes  in  the  media 
properties.  This  is  reasonable  because  if  the  temporal  distribution  of  events  is  small,  we  can 
assume  that  the  properties  of  the  media  have  not  changed  significantly,  and  if  the  spatial 
distribution  is  also  small  then  the  changes  in  the  travel  times  introduced  by  variations  in 
velocity  structure  is  minimized  as  well  (Fremont  and  Malone,  1987). 

The  second  possible  cause  of  waveform  variation  that  can  be  determined,  if  not  controlled, 
is  the  focal  mechanisms  of  the  events.  Fremont  and  Malone  (1987)  suggest  that  variations 
in  focal  mechanisms  usually  affects  the  waveform  shape  more  than  the  arrival  times  of  the 
P  and  S  waves  and  therefore,  do  not  constitute  as  strong  a  constraint  as  those  imposed  by 
the  spatial  and  temporal  correlations.  In  our  study,  we  will  show  that  the  focal  mechanisms 
for  the  Larderello  cluster  exhibit  a  high  degree  of  similarity  and  therefore,  no  degradation 
of  the  waveform  is  expected  from  this  variation  source. 

If  the  events  in  the  cluster  exhibit  similar  waveforms,  then  we  can  use  cross-correlation 
techniques  to  determine  the  degree  of  similarity  between  them,  and  then  either  eliminate 
outliers  completely  or  weight  them  appropriately  in  the  relative  location  inversion.  By 
identifying  and  controlling  variations  from  other  sources,  we  maximize  the  likelihood  that 
any  remaining  variations  are  attributable  only  to  changes  in  hypocenter  coordinates. 
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Interest  in  the  occurrence  of  events  showing  identical  or  nearly  identical  waveforms  has 
increased  in  recent  years  as  high  precision  location  techniques  are  developed.  Multiple 
events  with  a  strong  similarity  in  waveform  have  been  classified  into  ’’families”  or  termed 
doublets,  triplets  or  multiplets  according  to  the  number  of  similar  events  observed  (Geller 
and  Mueller,  1980;  Frankel,  1982;  Poupinet  et  al. ,  1984;  Fremont  and  Malone,  1987).  An 
earthquake  cluster  exhibiting  these  characteristics  is  an  excellent  candidate  for  relocation 
using  differential  travel  time  techniques. 


The  Relocation  Procedure 

The  location  of  source  parameters  by  the  method  of  differential  travel  times  uses  least 
squares  inversion  to  iteratively  solve  the  system  of  equations  £,•  -  F,-  =  e,-  ,  where  t;  is  the 
observed  travel  time  to  the  ith  station,  F,  is  the  calculated  travel  time  to  the  ith  station,  and 
et-  is  the  observational  error  in  t,  (Flinn,  1965;  Jordan  and  Sverdrup,  1981;  Menke,  1989). 

The  multiple  event  relocation  method  used  here  solves  a  non-linear  system  which  is  of  the 
form: 


d  =  A(h)  +  Bt  +  n  (0.1) 

where;  the  data  vector  d  contains  the  absolute  and  differential  arrival  times  for  a  set  of 
events,  stations  and  phases,  the  model  vector  h  contains  the  hypocenter  coordinates  of 
the  events,  the  t  vector  contains  the  origin  times  of  the  events,  and  the  vector  n  contains 
all  errors  attributable  to  observation  irregularities  and  approximations  in  the  estimate  of 
the  seismic  velocity  model.  ”A”  is  a  nonlinear  operator  obtained  from  the  function  that 
determines  the  travel  time  of  a  seismic  ray  (based  upon  an  estimate  of  the  velocity  structure, 
the  estimated  hypocenter  locations  and  the  location  of  the  recording  stations) .  B  is  a  linear 
operator  derived  from  the  differential  origin  times  of  a  slave  event  relative  to  the  origin  time 
of  a  master  event  .  A  and  B  are  implicitly  defined  by  models  of  the  absolute  or  differential 
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arrival  times. 


For  absolute  arrival  times  our  model  is: 

Tesw  =  te  +  F(xe,  x„  Uw)  +  nfssw  (0.2) 

and  for  differential  arrival  time  referred  to  one  master  event  we  use 

STresw  =  te  +  F{xe ,  x„  Uw)  +  nf*w  (0.3) 

where  the  following  definitions  apply: 

•  subscripts  e,  s  and  w  refer  to  event,  station  and  wave  type,  respectively 

•  STresw  refers  to  differential  travel  times  for  event  e  and  wave  w  at  station  s  referred 
to  a  reference  event  r 

•  te  and  tr  refer  to  the  origin  times  of  a  particular  event  and  reference  event,  respectively 

•  X{  is  a  vector  containing  the  hypocenter  coordinates 

•  Uw  is  the  spatial  dependence  of  the  seismic  slowness  (reciprocal  of  velocity) 

•  F(xe ,  xs,  Uw)  is  the  function  defining  the  travel  time  of  a  seismic  ray  from  the  hypocen- 
tral  coordinates  of  event  e  to  the  station  s  as  determined  from  the  seismic  slowness 
function  (a  function  of  position) 

•  nl\sw  is  the  observational  error  as  it  pertains  to  event  e,  station  s,  and  wave  w . 

An  additional  source  of  error  exists  in  that  the  slowness  function  is  approximated  by  a 
one-dimensional  model.  This  approximation  affects  the  absolute  relocations,  but  has  been 
shown  to  have  limited  effects  on  the  relative  relocation  of  events  whose  true  locations  are 
separated  by  less  than  one  wavelength  (Rodi  et  al.,  1993). 

Estimates  of  the  source  parameters  (the  h  and  t  vectors  in  equation  (0.1)  are  found  by 
minimizing  the  misfit  between  the  observed  and  calculated  travel  times  by  using  the  Polak- 
Ribiere  variant  of  the  conjugate  gradient  method  (Rodi  et  al.,  1993). 
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The  differential  times  are  determined  by  calculating  the  degree  of  similarity  between  a 
master  and  slave  event  by  cross-covariance  and  equivalently,  cross-correlation  techniques. 
Given  two  time  series  x  and  y,  the  cross-covariance  function  is  defined  as: 


Jxy 


rt2-ti 

(*)=  / 

JU 


x(t')y(tf  -  t)dt 


(0.4) 


and  the  cross-correlation  function  (the  normalized  form  of  the  cross-covariance  function)  is: 


CXy(t)  — 


\J  SXx{Q)Syy  (0) 


(0.5) 


If  the  two  events  to  be  correlated  are  identical,  equation  (0.5)  becomes  the  auto-  covariance 
function.  The  differential  times  are  determined  by  measuring  the  difference  between  the 
maximum  peak  of  the  cross-covariance  function  and  the  maximum  peak  of  the  auto-  covari¬ 
ance  function.  The  resulting  differential  times  give  the  errors  in  the  picks  of  the  arrival  for 
the  slave  events  relative  to  the  master  events.  The  differential  travel  times  were  determined 
directly  from  the  cross-covariances,  where  as  the  correlation  coefficients  were  determined 
from  the  cross-correlation  function,  as  defined. 


0.2.2  The  Method  used  to  Determine  Microearthquake  Focal  Mecha¬ 
nisms 

Focal  mechanisms  for  the  Larderello  microearthquakes  were  determined  using  a  computer 
algorithm  developed  by  Stewart  Guinn  and  L.T.  Long  (1977).  This  algorithm  uses  forward 
modelling  to  predict  the  polarity,  or  ’’sense”  of  the  first  motion  P  wave  arrivals  for  a  set 
of  seismic  stations.  The  first  motions  are  calculated  based  upon  P  wave  radiation  patterns 
from  a  dislocation  source  on  a  fault  plane  with  known  strike,  dip  and  rake.  Valid  solutions 
include  all  fault  planes  that  produce  a  match  between  the  calculated  first  motions  and  the 
observed  ones  (Hermann,  1975;  Guinn  and  Long,  1977). 
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The  fault  planes  are  determined  by  conducting  a  systematic  search  of  all  possible  orienta¬ 
tions  of  the  local  stress  ellipse,  which  is  related  to  the  orientation  of  the  fault  plane  through 
its  normal  and  slip  vectors  (Aki  and  Richards,  1980).  The  search  begins  by  orienting  the 
null  stress  axis  (the  B  axis)  in  the  vertical  direction  and  then  rotating  the  compressional 
and  tensional  axes  (P  and  T,  respectively)  around  the  B  axes  at  a  preset  angle  of  rotation. 
The  algorithm  then  calculates  the  strike  and  dip  of  the  two  nodal  planes  and  two  rakes  as¬ 
sociated  with  each  orientation  of  the  stress  ellipse,  and  based  on  these,  predicts  the  ’’sense” 
of  the  P  wave  first  motions  at  each  station.  In  order  to  accommodate  the  potential  for 
incorrect  first  motion  determinations,  the  user  may  specify  the  conditions  under  which  a 
solution  is  accepted  by  allowing  a  set  number  of  inconsistent  first  motions  in  the  solution 
set.  The  final  output  consists  of  the  azimuth  and  plunge  of  each  axis  and  the  strike,  dip 
and  rake  for  each  nodal  plane  that  produces  a  match  (Guinn  and  Long,  1977). 

All  solutions  which  match  the  observed  data  are  recorded  as  equally  valid  solutions  with 
three  ’’axis  domains”  defined.  For  well  constrained  problems,  the  domains  are  small.  For 
example,  if  50  solutions  are  found  that  match  the  data,  a  well  constrained  solution  set  will 
have  all  50  B  axes  located  within  a  small  area,  and  likewise  for  the  other  axes.  If  the  solution 
set  is  poorly-  constrained  then  the  axis  domains  will  define  large  and  perhaps  convoluted 
shapes  (Guinn  and  Long,  1977). 

0.2.3  Results 

The  1993  Larderello  microearthquake  cluster  is  an  appropriate  choice  for  the  application 
of  multiple  event  relocation  using  differential  travel  times  because  it  meets  the  criteria  dis¬ 
cussed  earlier.  The  cluster  occurred  in  a  spatially  restricted  zone  near  the  northern  boundary 
of  the  Larderello  geothermal  region.  Although  the  initial  hypocentral  locations  suggested 
diffuse  seismicity  over  a  volume  of  approximately  5  cubic  kilometers,  an  evaluation  of  the 
waveforms  indicated  a  high  degree  of  similarity,  suggesting  that  the  true  locations  are  more 
closely  located  than  that  suggested  by  the  single  event  inversion.  Temporally,  the  cluster 
occurred  over  the  course  of  a  few  days,  with  origin  time  intervals  ranging  from  minutes  to 
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hours  (See  Figure  0-3).  The  distribution  of  magnitudes,  in  combination  with  waveform 
similarity,  suggests  that  the  events  can  be  classified  into  multiplets,  which  supports  the  use 
of  cross-correlation  techniques. 

The  cluster  was  recorded  by  the  Larderello  seismic  network  which  is  a  local  array  of  thirty- 
one  short  period  seismometers.  ENEL  provided  MIT  with  seismic  records  of  200  events 
recorded  by  twenty-one  of  the  thirty-one  stations.  Nineteen  of  these  stations  were  equipped 
with  seismometers  that  recorded  the  vertical  component  of  ground  motion  only.  The  re¬ 
maining  two  stations  recorded  all  three  components,  i.e.  north,  east  and  vertical. 

The  events  included  in  this  study  were  chosen  such  that  there  was  a  high  percentage  of 
stations  actively  recording  during  that  time  period.  In  addition,  100  consecutive  events 
were  chosen  rather  than  100  randomly  distributed  ones.  This  was  done  so  that  if  it  was 
later  discovered  that  a  high  degree  of  correlation  existed  between  the  temporal  and  spatial 
distribution  of  events,  then  the  variation  in  velocity  structure  between  consecutive  events 
would  be  as  small  as  possible,  improving  the  likelihood  that  the  waveform  correlation  would 
be  high  as  well. 

Relocation  of  the  Larderello  cluster  was  accomplished  in  three  phases  beginning  with  the 
determination  of  the  absolute  and  differential  times,  followed  by  inversion  for  event  source 
parameters  (i.e.  hypocentral  locations  and  origin  times),  and  finally  calculation  of  focal 
mechanisms.  The  following  discusion  outlines  the  methods  used,  detail  the  results,  and 
present  an  interpretation  in  light  of  other  geologic  and  geophysical  data. 

Determining  differential  travel  times  for  the  Larderello  earthquakes  was  completed  in  seven 
stages: 

•  P  and  S  phase  picking 

•  Identifying  multiplets 

•  Determining  P  and  S  window  lengths  and  filtering  parameters 
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•  Choosing  master  events 


•  Choosing  stations  to  be  included  in  the  inversion 

•  Performing  the  cross-correlations 


P  and  S  arrival  times  were  identified  and  marked  in  the  headers  of  86  events  recorded  at 
21  stations.  The  sense  of  the  P  wave  first  arrival  was  also  recorded,  in  anticipation  of  the 
determination  of  focal  mechanisms.  Figure  0-4  presents  records  from  one  station.  The 
sampling  interval  for  all  records  was  0.004  seconds.  Preliminary  estimates  of  0.02  and  0.04 
seconds  were  assigned  for  the  standard  deviation  of  the  P  and  S  arrival  time  picks  (referred 
to  as  ’absolute’  arrival  times). 

During  the  P  and  S  phase  picking,  two  characteristic  waveforms  were  identified.  These 
waveforms  were  designated  as  Types  I  and  II  and  appear  to  be  characteristic  of  two  multi- 
plets  occurring  within  the  larger  cluster.  Referring  to  Figure  0-4,  note  that  events  147  and 
148  are  similar  and  of  Type  II  and  Events  149  and  150  are  similar  and  of  Type  I.  Event  184 
is  also  of  Type  II,  but  the  similarity  is  less  pronounced  than  for  the  others.  (The  numerical 
designation  indicates  the  relative  order  in  origin  times  of  the  events).  The  waveforms  for 
events  147  and  148  are  similar  in  their  full  length,  as  are  the  records  from  events  149  and 
150.  The  magnitudes  of  these  events  are  not  identical  in  the  pairs,  but  are  similar  (e.g. 
1.14  vrs.  1.34  for  events  149  and  150,  respectively).  Despite  variations  in  magnitude,  the 
similarities  in  waveform  persist,  and  we  note  that  in  general,  the  events  which  could  be 
classified  as  either  Type  I  or  II  were  found  to  have  a  magnitude  between  0.76  and  1.55. 
Figures  0-5  and  0-6  show  the  cross-covariances  of  these  event  records  with  an  event  of  each 
type  (II  and  I,  respectively).  These  cluster  characteristics  are  similar  to  those  found  by 
other  research  groups  investigating  multiplets  from  other  earthquake  clusters  (Fremont  and 
Malone,  1987;  Pechmann  and  Thorbjarnardottir,  1990;  Deichmann  and  Garcia-Fernandez, 
1992). 

Window  lengths  for  the  P  and  S  phases  were  determined  for  each  station  by  estimating 
the  maximum  possible  time  duration  that  would  exclude  any  portion  of  a  later  arrival.  For 
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example,  the  window  length  for  Type  II  events,  such  as  147  (Figure  0-4)  needed  to  have  an 
offset  less  than  0.35  seconds  after  the  P  arrival  in  order  to  exclude  the  second  arrival.  The 
full  length  of  the  windows  are  a  compromise  between  the  need  to  exclude  these  later  arrivals 
and  the  need  to  include  enough  cycles  of  the  desired  arrival  to  make  the  correlations  valid. 
In  addition,  the  actual  window  lengths  for  each  station  and  phase  were  determined  based 
upon  an  evaluation  of  the  amplitude  spectra  for  events  recorded  at  the  station.  In  all  cases, 
the  window  lengths  were  between  0.5  and  0.8  seconds  long  depending  upon  the  station. 
Figure  0-7  shows  an  example  of  an  event’s  amplitude  spectrum  at  station  CORN.  The 
majority  of  spectral  energy  falls  between  5  and  15  Hertz  with  a  marked  increase  between  8 
and  13  Hertz.  Window  lengths  of  0.5  seconds  allow  four  cycles  of  an  8  Hertz  wave.  Longer 
window  lengths  were  used  whenever  the  records  for  a  station  indicated  the  absence  of  a 
second  arrival  within  that  interval,  so  that  the  lower  frequencies  would  be  represented  in 
the  correlation.  Prior  to  performing  the  correlations,  the  full  waveforms  were  bandpass 
filtered  (corners  of  5  and  20  Hertz)  and  a  cosine  taper  was  applied  to  the  window. 

Master  events  were  selected  on  the  basis  of  clarity  and  impulsiveness  of  P  and  S  arrivals, 
low  noise  and  good  coverage  by  stations  (i.e.  masters  recorded  at  as  many  stations  as 
possible).  Several  attempts  were  made  to  use  only  two  masters  for  the  entire  cluster  of 
events,  however,  an  unacceptable  proportion  of  low  correlations  required  the  use  of  multiple 
masters.  So,  the  cluster  was  divided  into  four  subsets,  each  with  2  masters,  one  of  each 
type.  For  example,  events  150  and  156  were  used  as  masters  for  events  141  through  156. 
For  events  116  through  140,  a  strong  type  II  master  was  not  identified,  thus,  event  126  was 
included  as  a  secondary  master  for  that  group.  Figures  0-8  and  0-9  present  the  master 
events  at  station  POMA. 

Stations  used  in  the  inversion  exhibited  the  following  characteristics;  most  of  the  master 
events  and  many  of  the  regular  events  were  recorded,  and  the  signal-to-noise  ratio  was 
acceptable.  Additionally,  distance  and  azimuthal  coverage  were  also  used  as  criteria.  The 
stations  that  were  excluded  exhibited  poor  signal-to-noise  ratios  or  recorded  only  a  limited 
number  of  events.  The  final  stations  chosen  were:  CLSV,  CORN,  FRAS,  MINI,  MONV, 
POMA  and  SCAP  (See  Figure  0-2). 


14 


November  20,  1997  19:01 


Differential  travel  times  for  each  phase  (P  and  S)  were  determined  for  each  of  the  86  events 
at  each  of  the  stations.  The  following  master-slave  event  pairs  were  used;  events  126  and 
131  were  masters  for  events  116-140,  events  150  and  156  were  masters  for  events  141-156, 
events  171  and  186  were  masters  for  events  157-187  and  events  189  and  191  were  masters 
for  events  188-218. 

Differential  travel  times  were  calculated  for  each  phase  using  the  following  equation: 

dt  =  ttmaster  tts[ave  C picking  t waveform  (0.6) 


where, 


master  —  master  TO 


master 


t^slave  —  dtslave  T0slave 


In  practice,  the  use  of  the  origin  times  is  avoided  by  introducing  a  common  reference  time 
for  both  events  and  normalizing  the  arrival  times  to  the  master  events  arrival  times. 

Figure  0-11  shows  a  schematic  diagram  of  the  elements  of  the  differential  time.  For  ease 
of  illustration,  we  choose  to  demonstrate  the  process  of  determining  differential  times  using 
a  master  and  slave  pair  that  are  coincident  in  time  (as  indicated  in  the  figure  by  their 
coincident  origin  times,  TO).  The  master  event  arrival  times  show  a  unique  moveout  curve 
for  any  particular  array  of  stations.  (The  array  we  are  using  here  is  a  linear  array  which 
we  choose  for  convenience  in  order  to  show  a  simply  shaped  moveout  curve).  The  slave 
event  arrival  time  moveout  curve  exhibits  a  slightly  different  shape,  which  is  an  expres¬ 
sion  of  the  difference  in  location  between  the  master  and  the  slave  events.  When  the  two 
events  are  identical  in  all  ways  except  location  and  we  declare  (for  this  example)  that  the 
picking  error  is  zero,  then  the  epsilons  in  equation  (0.6)  are  zero  and  we  can  derive  the 
differential  times  directly  from  the  absolute  time  picks.  Usually,  however,  picking  error  is 
not  negligible  and  there  are  some  variations  in  the  waveform  due  to  the  change  in  location. 
We  assume  that  there  are  no  variations  in  focal  mechanism  and  source  time  functions  and 
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so  all  waveform  variations  are  due  to  changes  in  hypocenter  location  only.  If  the  events 
are  sufficiently  close  together  then  the  waveform  variations  will  be  small  and  we  can  use 
correlation  techniques  to  remove  any  temporal  offsets  (epsilons)  caused  by  the  slight  wave¬ 
form  differences.  Additionally,  the  correlation  will  also  remove  any  picking  inconsistencies 
by  using  the  pick  on  the  master  event  (per  station)  as  a  reference.  It  will  not,  however, 
remove  picking  inconsistencies  between  stations. 

The  temporal  offsets  (also  refered  to  as  correlation  offsets)  were  determined  by  cross¬ 
covariance  of  a  P  or  S  window  from  a  slave  event  and  a  corresponding  P  and  S  window 
from  a  master  event.  Figures  0-10  through  1  0-14  illustrate  the  determination  of  the  off¬ 
sets  from  the  cross-covariances.  The  time  lag  between  the  cross-covariance  peak  and  the 
auto-covariance  peak  is  the  correlation  offset. 

The  correlation  offsets,  taken  as  a  set  of  data,  are  expected  to  exhibit  the  statistical  char¬ 
acteristics  of  a  gaussian  distribution  with  zero  mean  and  variance  (Jordan  and  Sverdrup, 
1981;  Rodi  et  al.,  1993)  Figures  0-15  and  0-16  present  examples  of  the  distributions  of 
cross-  correlation  offsets  for  P  and  S  waves  for  the  86  events  from  the  Larderello  cluster 
(as  recorded  at  station  POMA).  The  distribution  of  the  P  wave  offsets  clearly  shows  the 
expected  characteristics,  while  the  distribution  of  the  S  wave  offsets  suggest  that  the  orig¬ 
inal  S  wave  arrival  time  picks  were  affected  by  a  non-random  influence.  This  influence  is 
most  likely  related  to  the  masking  of  the  S  wave  first  arrivals  by  P  wave  coda  and  P-to-S  or 
S-to-P  conversions.  For  some  of  the  stations  and  events,  earlier  arrivals  were  superimposed 
on  the  S  wave  first  arrivals  making  them  difficult  to  pick.  For  those  stations  and  events, 
the  expected  S  wave  arrival  times  were  anticipated  based  upon  a  set  of  clear  arrivals  and 
three-component  records  (when  available).  These  reference  S  waves  were  associated  with 
a  particular  peak  in  the  waveform  and  then  traced  through  the  rest  of  the  records.  This 
technique  would  maximize  the  likelihood  that  the  true  S  arrival  would  occur  within  the 
chosen  window,  allowing  the  cross-correlation  to  detect  it,  but  the  bias  would  then  appear 
as  a  non-random  effect  in  the  distribution  of  the  correlation  offsets. 

The  standard  deviations  of  the  P  and  S  differential  travel  times  were  estimated  by  taking,  as 
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a  baseline,  values  of  0.004  seconds  for  the  P  waves  and  0.008  seconds  for  the  S  waves.  These 
values  were  then  modified  for  each  master-slave  pair  by  dividing  by  the  cross-correlation 
coefficient,  so  that  pairs  showing  the  highest  degree  of  correlation  would  have  the  smallest 
possible  standard  deviations. 
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Depth  Below 
sea  level 

P  Wave 
velocity  (m/s) 

S  Wave 
velocity  (m/s) 

0 

3,500 

2,000 

300 

5,500 

3,200 

3,000 

6,000 

3,500 

Table  0.1:  1-D  velocity  structure  model  (version  0). 

Hypocenter  Inversion 

The  multiple  relative  location  inversion  involved  two  steps: 

•  Developing  starting  models  for  the  hypocenters  and  1-D  velocity  structures 

•  Inverting  for  the  full  set  of  locations,  evaluating  the  results  and  estimating  errors 

•  Evaluating  errors 

The  starting  model  for  the  hypocenter  locations  was  based  upon  ENEL’s  initial  hypocenter 
estimates  from  the  single  event  inversion.  All  events  were  assigned  the  same  guess  of  25,000 
meters  north  of  the  local  coordinate  system  origin  and  33,000  meters  east  of  the  origin. 
This  starting  model  would  eventually  be  modified  several  times  to  reflect  the  results  from 
subsequent  inversion  trials.  Several  different  velocity  models  were  used  in  the  inversion. 
These  models  evolved  as  additional  information  became  available  regarding  known  geologi¬ 
cal  characteristics  of  the  region.  The  development  and  evolution  of  these  models  was  based 
on  available  geological  information  which  will  be  discussed  in  detail  in  a  later  section  of  this 
paper  (Tables  (0.1),  (0.2),  and  (0.3). 

The  hypocenters  and  origin  times  were  determined  using  the  non-linear,  conjugate  gradient 
inversion  algorithm  described  earlier.  Various  starting  models  were  tried  until  the  best  fit 
to  the  observed  data  were  obtained.  The  final  results  were  found  using  a  starting  model  of 
25,000  m  north,  35,000  m  east  and  4600  m  depth  and  version  2  of  the  velocity  model  (0.3). 
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Depth  Below 
sea  level 

P  Wave 
velocity  (m/s) 

S  Wave 
velocity  (m/s) 

0 

3,500 

2,000 

500 

4,800 

2,700 

1,000 

5,900 

3,200 

4,800 

4,900 

2,800 

Table  0.2:  1-D  velocity  structure  model  (version  1). 


Depth  Below 
sea  level 

P  Wave 
velocity  (m/s) 

S  Wave 
velocity  (m/s) 

4,500 

2,360 

4,800 

2,800 

500 

5,900 

3,200 

6,100 

3,500 

4,900 

2,865 

4,200 

6,100 

3,500 

Table  0.3:  1-D  velocity  structure  model  (version  2). 

Locations  are  shown  in  Figures  0-17  through  0-20,  and  the  distribution  of  the  residuals 
for  the  relative  travel  times  are  presented  in  Figure  0-21. 

To  evaluate  the  results  in  terms  of  geometrical  patterns  we  must  consider  two  possible 
sources  of  bias;  1)  hypocenters  with  large  residuals,  and  2)  offsets  in  the  locations  of  a 
subset  relative  to  another  subset  due  to  the  use  of  multiple  masters. 

Scatter  introduced  by  poorly  fit  locations  was  eliminated  by  culling  outliers  on  the  basis  of 
residual  travel  times,  data  standard  deviation  estimates  and  error  ellipses  from  the  inversion. 
All  events  whose  travel  time  residuals  exceeded  three  times  the  standard  deviations  of  the 
data  were  excluded  from  the  geometry  analysis.  This  is  a  stricter  criteria  than  using  the 
standard  deviations  of  the  residuals  and  was  chosen  so  that  only  the  best  fit  data  was  used 
in  the  geometry  analysis. 
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Errors  due  to  multiple  masters  exist  in  the  set  of  locations  because  the  constraint  placed 
on  the  inversion  by  the  differential  travel  times  between  master-master  pairs  was  relatively 
weak.  Any  geometric  pattern  in  the  locations  may  therefore  be  obscured  by  errors  in  the 
relative  locations  between  subsets.  It  is  necessary  then,  to  eliminate  the  effects  of  these 
errors  by  separating  the  subsets  and  viewing  them  inpedendently.  For  each  subset,  one 
master  was  used  as  the  dominant  reference  event  because  it  was  observed  that  one  or  the 
other  of  the  masters  produced  higher  correlations.  Separating  on  the  basis  of  masters  ensures 
that  master  related  scatter  is  removed.  The  separated  and  culled  results  are  presented  in 
Figures  0-22  through  0-25. 

Two  of  the  subsets  (events  157-187  and  188-218)  show  distinct  lineations.  The  two-dimensional 
confidence  regions  (for  epicentral  locations)  and  confidence  intervals  (for  depth)  give  an  in¬ 
dication  of  the  likelihood  of  the  apparent  patterns.  Figures  0-26  through  0-29  present 
some  typical  estimated  confidence  regions  for  the  subsets  of  interest. 

The  absolute  locations  have  been  estimated  from  the  master  events,  for  example,  the  ab¬ 
solute  location  of  event  191  is  25,808  m  north,  35,033  m  east  and  5,070  m  in  depth.  The 
uncertainties  in  absolute  location  for  this  event  are  625  m  (semi-major  axis),  503  m  (semi¬ 
minor  axis),  and  1,049  m  in  depth. 

The  pattern  and  shape  of  the  error  ellipses,  as  determined  by  the  relative  event  locations 
for  subset  number  4  (events  188-218),  indicate  that  the  apparent  lineation  is  not  an  ar¬ 
tifact  of  uncertainties  in  the  inversion  results.  Event  157-187  exhibit  greater  uncertainty 
in  their  locations,  however,  the  common  overlap  of  the  elipses  defines  a  very  small  region. 

It  is  unlikely  that  the  northwest-southeast  trend  of  the  event  locations  can  be  reasonably 
attributed  to  this  uncertainty.  In  light  of  this,  the  orientation  of  these  two  lineations  should 
be  compared  with  the  the  known  geological  structures  and  any  other  possible  site  specific 
associations.  Patterns  in  relative  depth  distribution,  however,  are  not  well  resolved,  and 
therefore,  any  conclusions  regarding  depth  should  be  derived  from  the  absolute  locations. 
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Focal  Mechanisms 


Focal  mechanisms  were  determined  for  17  events  in  the  Larderello  microearthquake  cluster 
using  first  motion  P  wave  data.  The  clarity  of  the  P  wave  arrivals  was  excellent  at  most 
stations  and  no  problems  were  encountered  discerning  the  sense  of  the  motion,  as  evidenced 
in  Figures  0-8,  0-9,  and  0-10.  Several  stations  exhibited  variations  in  the  direction  of 
motion  indicating  that  rays  from  these  stations  were  emerging  near  the  nodal  planes.  These 
stations  provided  critical  information  and  additional  constraint  on  the  final  fault  plane 
solutions.  The  solutions  for  the  17  events  indicated  little  variation  between  mechanisms. 
Figure  0-30  represents  examples  of  the  extreme  cases  observed.  An  effort  was  made  to  use 
as  many  stations  as  possible  and  select  events  with  the  best  coverage  available,  however, 
the  azimuthal  and  distance  coverage  was  impacted  by  a  lack  of  stations  in  the  northwestern 
quadrant  and  the  location  of  the  swarm  near  the  northern  boundary  of  the  network.  Events 
for  which  focal  mechanisms  could  be  determined  were  recorded  by  between  13  and  16 
stations,  generally  including  CLSV,  CORN,  CRBE,  FRAS,  FROS,  LAGO,  LURI,  MBAM, 
MDSV,  MGUI,  MINI,  MONV,  POMA,  SCAP,  TRAV,  VALE,  and  MLUC  (Figure  0-17). 

The  number  of  valid  solutions  identified  for  each  of  the  events  varied  according  to  how  many 
stations  were  available  and  the  particular  set  of  first  motions,  however,  all  events  exhibited 
well  defined  domains  for  the  three  principle  axes  of  the  stress  tensor.  Focal  mechanisms  were 
calculated  using  a  search  increment  of  10  degrees,  as  recommended  by  Guinn  and  Long, 
and  the  number  of  inconsistent  first  motions  was  limited  to  one  or  zero  (Guinn  and  Long, 
1977).  Solutions  for  all  events  indicate  reverse  faulting  mechanisms  ranging  from  nearly 
pure  reverse  to  predominantly  reverse  with  some  strike-slip  component.  Average  solution 
parameters  for  two  representative  events  are  listed  in  Tables  (0.4)  and  (0.5). 
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■1 

P  Axis 

T  Axis 

B  Axis 

Event 

Azimuth 

Plunge 

Azimuth 

Plunge 

Azimuth 

Plunge 

200 

221 

86 

128 

51 

315 

40 

187 

61 

74 

196 

23 

325 

75 

Table  0.4:  Principle  axes  of  stress  for  events  200  and  187. 


Event 

Strike 

Dip 

Rake 

NE  Dipping 
Nodal  Plane 

200 

168 

66 

147 

187 

173 

32 

12 

SW  Dipping 
Nodal  Plane 

200 

272 

61 

28 

187 

319 

63 

73 

Table  0.5:  Strike,  dip  and  rake  of  nodal  plane  for  events  200  and  187. 

0.2.4  Interpretation 

Geologic  Setting 

The  Tuscany  region  of  Italy  is  associated  with  the  Tyrrhenian-Apennine  orogenic  system 
which  trends  northwest  through  the  Italian  peninsula  and  marks  the  location  of  a  sharply 
curved  thrust  belt  that  extends  south  and  then  west  through  Calabria  and  Sicily.  Structures 
within  this  region  reflect  at  least  three  separate  phases  of  tectonic  disturbances,  which  are 
recognized  by  remnant  features  in  both  the  sedimentary  formations  as  well  as  in  some  of 
the  basement  units.  The  primary  structures  of  the  Apennines  formed  during  the  Apennine 
Orogeny  (middle-Apline  Orogeny,  middle  Eocene  to  Oligocene,  30-50  ma)  (Carmignani  and 
Kligfield,  1990;  Cameli  et  a/.,  1993). 

The  first  disturbance,  the  Eocene-Oligocene  phase  of  deformations,  is  referred  to  by  Carmignani 
and  Kligfield  (1990)  as  the  D1  phase  and  is  represented  by  structural  features  created  by 
northeast-southwest  oriented  compression  during  the  initial  collision  of  the  Corsica-Sardinia 
microplate  with  the  Adriatic  microplate  .  The  D1  phase  created  nappe  structures  by  dis¬ 
placing  preexisting  sedimentary  formations  along  detachment  horizons  (Carmignani  and 
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Kligfield,  1990). 


The  sedimentary  formations  involved  in  the  nappe  building  are  identified  as  the  Ligurid 
complex  and  the  Tuscan  Nappe.  The  Ligurid  complex  is  composed  of  Jurassic  to  Paleogene 
sediments,  specifically  flysch,  pelagic  sediments,  radiolarites  and  ophiolites.  The  Tuscan 
nappe  are  composed  of  Triassic  to  Paleogene  flysch,  radiolarites,  marls,  limestones  and 
anhydrites.  These  two  units  are  overlain  by  neogene  sediments  (Cameli  et  al. ,  1993;  Batini 
et  al .,  1985a). 

The  allocthonous  Tuscan  nappe  is  discontinuous  due  to  extensional  thinning  and  in  some 
locations,  is  underlain  or  replaced  by  overthrust  wedges  of  the  autochthonous  Tuscan  series 
(referred  to  as  the  ’’tectonic  wedges”)  and  underlying  core  complex  metamorphic  formations 
(Carmignani  and  Kligfield,  1990;  Block,  1991). 

The  Tectonic  slices  is  part  of  the  Monteciano-Roccastrada  Unit  (or  Tuscan  Metamorphic 
sequence).  This  Mesozoic-Paleozoic  sequence  is  composed  of  three  units,  1)  quartz  metacon¬ 
glomerates,  quartzites  and  phyllites,  2)  phyllites  and  quartzites,  and  3)  micaschists.  Under¬ 
lying  the  Monteciano-Roccastrada  formation  lies  a  pre-Alpine  gneiss  complex  interbedded 
with  amphibolites  (Ordovician-Silurian  age).  Metamorphism  affecting  this  unit  originated 
in  the  Hercynian  Orogeny  (approximately  300  m.a.)  and  included  the  formation  of  NE- 
SW  trending  fractures  (Batini  et  al.,  1985a;  Cameli  et  al.,  1993;  Carmignani  and  Kligfield, 
1990).  The  Gneiss  and  Monteciano-Roccastrada  units  are  often  referred  to  in  the  literature 
as  being  representative  of  the  lower  plate  and  the  term  ’’basement”  has  been  applied  to 
both. 

The  compressional  stress  that  lead  to  the  displacement  of  the  Ligurid  and  Tuscan  sedimen¬ 
tary  sequences  onto  the  core  complex  was  followed  by  a  period  of  post  collisional  exten¬ 
sion  (referred  to  by  Carmignani  and  Kligfield  as  the  D2  deformation).  The  late  Tortonian 
(Miocene)  D2  tensional  stresses  created  boudinage  structures  within  the  cherty  limestones  of 
the  metamorphic  complex,  attesting  to  the  extensional  affects  at  mid-crustal  depths  (below 
3  km).  Effects  of  the  D2  deformation  on  the  upper  crust  can  be  seen  in  the  northwest- 
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erly  trending  extensional  normal  faults.  Vertical  displacement  along  these  fractures  created 
Horst-Graben  structures  into  which  sediments  deposited  during  the  Miocene  and  Pliocene 
(Neogene  sediments)  (Cameli  et  al .,  1993;  Carmignani  and  Kligfield,  1990). 

The  third  episode  of  tectonic  disturbance  includes  Pliocene  and  later  age  magmatism  and 
volcanics.  The  D2  extensional  deformation  creating  the  NW-SE  fractures  migrated  eastward 
across  the  Italian  peninsula  during  the  Pliocene  and  Quaternary  times.  The  existence  of  the 
faulting  provided  flow  pathways  for  later  magmatism  (Puxeddu,  1984).  This  magmatism  is 
thought  to  be  responsible  for  the  formation  of  the  Quaternary  volcanics  located  in  Tuscany 
as  well  as  numerous  granitic  intrusions  emplaced  throughout  the  geothermal  regions.  It  has 
been  suggested  that  most  of  Tuscany  is  underlain  by  granitic  type  intrusions.  Petrologic 
data  suggest  the  age  of  origin  of  the  accessible  formations  to  be  2. 5-3.5  ma.  (Batini  et  al ., 
1985a). 

The  velocity  model  used  in  the  relocation  was  developed  using  information  from  several 
sources.  P  and  S  wave  velocities,  and  their  ratios  (Vp,  Vs,  and  Vp/Vs)  were  based  on 
results  of  a  seismic  reflection  survey  conducted  by  Batini,  Duprat  and  Nicolich  (1985), 
and  on  joint  hypocenter-velocity  inversion  from  local  earthquake  data  conducted  by  Block 
(1991).  The  thicknesses  for  the  layers  were  based  on  the  seismic  reflection  results  and  on 
detailed  contour  maps  of  the  primary  units  (Batini  et  al.,  1985a, b;  ENEL-Unita  Nazionale 
Geotermica,  1988a, c,b;  Societa  Elaborazione  Cartographiche,  1978a, b,c,d). 

Batini,  Duprat  and  Nicolich  (1985)  identified  several  seismically  reflective  horizons  which 
correlate  well  with  the  Liguride,  Tuscan,  and  ” Basement”  units.  These  horizons,  identified 
as  the  ”L”,  ”T”  and  ”B”,  respectively,  were  associated  with  velocities  ranging  from  4,000 
m/s  to  6,100+  m/s,  with  Vp/Vs  ratios  ranging  from  1.7  to  1.9,  depending  upon  the  forma¬ 
tion  (Block,  1991;  Batini  et  al.,  1985b).  In  addition,  another  reflector,  known  as  the  ”K” 
horizon,  has  been  identified.  It  is  regionally  extensive,  covering  large  areas  in  southern  Tus¬ 
cany  and  varies  in  depth  from  as  shallow  as  2.8  km  in  the  southern  Larderello  geothermal 
field  to  greater  than  9  km  near  Montalcinello  (  20  km  east)  and  other  locations  (Cameli 
et  al.,  1993;  ENEL-Unita  Nazionale  Geotermica,  1988a). 
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The  ”K”  horizon  is  located  within  the  gneissic  basement,  and  although  several  possible  in¬ 
terpretations  regarding  its  nature  have  been  proposed,  its  exact  origin  has  not  been  clearly 
established.  One  interpretation  suggests  that  the  ”K”  horizon  represents  a  zone  of  frac- 
turation  and  the  presence  of  fluid  filled  fractures  (Batini  et  al.,  1985b),  while  others  have 
suggested  that  it  represents  a  rheological  boundary  between  the  brittle  upper  unit  and  a 
ductile  lower  unit  (Cameli  et  al.,  1993).  The  seismic  reflection  results  delineate  sub- vertical 
normal  faults  which  cross-  cut  all  basement  units  including  the  gneissic  formation  and  the 
”K”  horizon  (Batini  et  al.,  1985a).  Additionally,  other  sources  indicate  that  normal  listric 
faults  flatten  into  the  ”K”  horizon  and  conclude  that  the  brittle-ductile  boundary  is  kine¬ 
matically  active  (Cameli  et  al.,  1993). 

In  the  southern  portion  of  the  Larderello  geothermal  field,  the  K  horizon  splits  into  sev¬ 
eral  individual  reflectors  which  are  relatively  shallow.  A  well,  designated  San  Pompeo  2, 
penetrated  through  the  Monteciano-Roccastrada  Unit  and  the  gneiss  unit,  into  the  upper 
portions  of  the  ”K”  horizon.  The  basement  units  were  found  to  be  highly  fractured  to 
a  depth  of  2,300  m,  where  an  impermeable  phyllite  and  micaschist  unit  was  encountered. 
At  2,930  m  drilling  progress  stopped  when  the  drill  string  encountered  a  zone  of  highly 
fractured  rock  with  hydrothermal  and  contact  metamorphism  mineral  assemblages,  fluid 
pressures  in  excess  of  240  bars  (24  MPa)  and  temperatures  in  excess  of  400C  (Cappetti 
et  al.,  1985). 


The  Larderello  Geothermal  Field 

The  Larderello  geothermal  field  is  a  vapor  dominated,  hydrothermal  type  field  producing 
steam  primarily  from  zones  associated  with  the  Triassic  anhydrites  at  the  base  of  the  Tuscan 
nappe  (Puxeddu  et  al.,  1977).  The  geothermal  source  at  Larderello  is  thought  to  be  a 
granitic  type  intrusion  at  depths  ranging  from  6  km  to  greater  than  25  km  (Batini  et  al., 
1985a;  Block,  1991).  In  support  of  this  theory  is  the  identification  of  aplitic  veinlets  in 
samples  from  deep  wells  (Minissale,  1991)  and  a  broad  seismic  velocity  low  extending  from 
depths  of  approximately  7  km  to  25  km,  approximately  the  depth  of  the  Moho  (Block,  1991; 


25 


November  20,  1997  19:01 


Carmignani  and  Kligfield,  1990). 


Cameli  et  al.  (1993)  report  that  the  occurence  of  seismicity  within  the  boundaries  of  the 
Larderello  geothermal  area  is  concentrated  between  1  and  7  km  in  depth,  with  a  peak 
occurance  between  3  and  5  km.  This  depth  range  correlates  well  with  the  ”K”  horizon, 
and  with  estimates  from  the  relocations  presented  here.  However,  any  interpretation  of  the 
cause  of  the  events  must  account,  not  only  for  the  horizontal  and  vertical  placement,  but 
also  for  the  focal  mechanism  parameters.  We  offer  an  evaluation  of  each  of  these  elements 
in  the  following  discussion. 

Only  two  previously  identified  geological  structures  are  discussed  in  the  literature  for  this 
area  and  depth  range;  the  ”K”  horizon  and  the  NW-SE  trending,  high  angle,  normal  faults. 
A  detailed  map  of  the  northern  boundary  of  the  Larderello  field  indicates  that  the  ”K” 
horizon  is  located  between  4.4  and  5  km  in  depth  and  that  the  area  is  cross-cut  by  two  of 
the  normal  faults  (ENEL-Unita  Nazionale  Geotermica,  1988a).  The  horizontal  location  of 
the  March  1993  cluster  in  relation  to  these  features  is  shown  in  Figures  0-31  and  0-32. 

A  review  of  the  focal  mechanisms,  spatial  geometry  and  temporal  distributions  suggest 
that  the  microearthquakes  in  this  cluster  differ  in  overall  character  from  those  studied 
recently  in  regions  of  active  volcamsm,  such  as  the  Mid- Atlantic  Ridge  (MAR)  and  at 
convergent  plate  boundary  volcanoes.  In  recent  literature,  the  study  of  naturally  occurring 
microearthquakes  has  been  limited,  in  general,  to  events  with  magnitudes  greater  than  3.0 
due  to  limitations  imposed  by  the  detection  limits  of  the  seismic  networks.  However,  in  these 
regions  the  events  studied  have  been  attributed  to  several  specific  mechanisms,  including, 
dike  emplacement,  crustal  extension,  crustal  loading,  cooling  and  crystallization  of  magma, 
and  magma  withdrawal. 

In  the  case  of  dike  emplacement,  the  earthquakes  have  been  observed  to  track  the  advancing 
magma  front.  Relocation  of  the  Larderello  earthquakes  indicates  no  such  temporal  progres¬ 
sions  and  based  on  this,  and  on  the  lack  of  evidence  supporting  recent  or  active  migration  of 
magma  at  mid-crustal  depths,  we  consider  this  mechanism  as  being  an  unlikely  candidate 
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for  these  events. 


In  an  evaluation  of  earthquake  swarms  on  the  MAR,  Bergman  and  Soloman  (1990)  con¬ 
cluded  that  much  of  the  seismicity  detected  exhibited  normal  faulting  and  is  attributable 
to  crustal  extension,  but  that  this  conclusion  is  based  primarily  on  events  with  magnitudes 
above  the  detection  limit  of  the  global  seismic  network  (Mb=4.5).  Dzurisin  et  al  (1991) 
concluded  that  normal  faulting  and  subsidence  near  the  Medicine  Lake  Volcano,  in  North¬ 
ern  California  is  most  likely  due  to  crustal  loading  by  the  volcano  and  crustal  thinning  due 
to  basin  and  range  extension.  In  this  study  as  well,  the  magnitudes  of  the  events  ranged 
between  3.0  and  4.6,  and  the  focal  mechanisms  were  different  from  those  of  the  Larderello 
events. 

The  characteristics  of  the  Larderello  events  are  inconsistent  with  those  observed  at  other 
volcanic  areas,  however,  several  other  possible  alternatives  exist  which  take  into  account 
the  unique  nature  of  a  cooling  igneaous  body.  To  evaluate  these  factors,  we  must  clearly 
define  the  characteristics  which  must  be  accounted  for.  The  focal  mechanisms  are  the 
most  consistent  characteristic  of  the  events  and  therefore  we  consider  the  two  possible  valid 
solutions  identified  earlier.  The  first  is  that  the  events  occurred  on  one  or  more  northwest 
striking  faults  which  dip  approximately  60  degrees  to  the  northeast.  Rakes  ranging  from  28 
to  73  degrees  indicate  primarily  reverse  faulting  with  varying  degrees  of  left  lateral  strike 
slip.  The  second  alternative  is  that  the  events  occur  on  south-southeast  striking  faults  with 
varying  degrees  of  dip  to  the  west,  and  varying  degrees  of  reverse  and  strike-slip  components. 

The  horizontal  location  of  the  cluster  and  the  first  focal  mechanism  alternative  strongly 
suggests  an  association  between  the  events  and  the  northern-most  high  angle  fault,  and 
taken  along  with  the  depth  estimates,  also  indicates  a  possible  interaction  between  the 
faults  and  the  ”K”  horizon.  The  second  alternative  clearly  excludes  the  northern-most  high 
angle  fault  (due  to  an  inconsistent  dip  direction) ,  however,  a  potential  association  with  the 
”K”  horizon  still  exists. 

Several  pieces  of  evidence  favor  the  first  alternative.  First,  the  correlation  between  strike 
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and  dip  of  the  fault  plane  solution  and  that  of  the  existing  faults  in  the  gneissic  unit 
is  compelling,  however,  it  remains  to  be  explained  how  the  D2  normal  faults  might  be 
reactivated  in  reverse  faulting.  One  possibility  is  that  subsidence  or  downwarping  of  the 
“K”  horizon  is  occurring  due  to  cooling  and  or  a  drop  in  pore  pressure  due  to  steam 
migration.  Another  possibility  is  that  microearthquakes  represent  only  small  displacements 
and  therefore,  the  focal  mechanisms  for  these  events  are  most  reasonably  associated  with 
very  local  phenomena  rather  than  regional  ones.  The  seeming  inconsistency  between  the 
reverse  faulting  mechanism  and  the  regional  extension  is  therefore  not  a  critical  barrier  to 
this  alternative.  In  light  of  this,  we  might  consider  the  events  to  represent  displacement 
of  small  blocks  within  or  near  the  fault  zone  with  parallel  or  sub-parallel  surfaces.  Figures 
0-24  through  0-27  seem  to  support  a  connection  with  the  two  fault  directions  apparent  in 
the  gneissic  unit  as  observed  by  Carmignani  and  Kligfield  (1990). 

0.2.5  Discussion 

The  results  of  the  relocation  of  the  March  1993  Larderello  microearthquake  cluster  place  the 
events  near  the  northern  boundary  of  the  vapor-dominated  geothermal  reservoir.  Proximity 
of  the  cluster  to  a  Tortonian  age  normal  fault  and  the  seismic  reflector  referred  to  as  the 
”K”  horizon  suggest  a  complex  interaction  between  the  structures,  temerature  changes, 
and  the  hydrologic  conditions  at  the  geothermal  field.  I  propose  that  these  events  may  be 
associated  with  thermal  contraction  of  the  granitic  body.  I  make  this  proposal  based  on  the 
focal  mechanisms  indicating  reverse  faulting  and  the  spatial  association  with  pre-existing 
structures. 
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0.3  CONCLUSIONS 


Our  work  on  the  1993  microearthquake  swarm  at  the  Larderello  geothermal  field  suggests 
that  the  most  significant  obstacle  to  the  practical  application  of  correlation  techniques  such 
as  EGFs  or  relative  event  location  for  seismic  events  recorded  at  regional  distances  is  due 
to  the  degree  of  waveform  complexity  that  is  related  solely  to  scattering  affects  along  the 
path.  In  this  study  we  showed  that  dispite  similar  focal  mechanisms  and  a  regional  scale 
that  is  on  the  order  of  1  to  2  wavelengths,  the  complexity  of  the  waveforms  was  severe, 
particularly  for  S  waves.  We  showed  that  the  event  waveform  types  exhibited  at  least  two 
distinctly  different  characters,  thus  forcing  the  use  of  more  than  one  master  event. 

Our  initial  analysis  of  the  cluster  size  suggested  that  for  P  waves  the  pertinent  wavelengths 
were  below  1220  meters,  based  upon  a  P  wave  velocity  of  6100  m/s,  with  a  mean  wavelength 
of  approximately  700  meters.  The  pertinent  wavelength  scale  for  S  waves  was  below  700 
meters  so  we  expected  to  have  more  difficulty  with  the  S  wave  correlation  coefficients. 
Although  the  size  of  the  overall  cluster  exceeded  these  bounds,  the  smaller  subsets  were  of 
the  required  size  (that  being  approximately  one  wavelength). 

Another  difficulty  arose  with  the  application  of  windowing  around  an  absolute  arrival  time 
pick.  This  method  was  applied  in  an  effort  to  bracket  the  true  arrival  time  and  to  reduce 
the  computing  operations  by  limiting  the  lengths  of  the  data  sequences  being  correlated. 
The  existence  of  strong  scatterers  introduced  a  significant  difficulty  in  choosing  the  correct 
arrival  and  limited  the  window  size,  thus  potentially  affecting  the  accuracy  of  the  correlation 
coefficient. 

The  difficulties  discussed  above  were  partially  mitigated  by  a  very  careful  and  detailed 
analysis  of  correlation  coefficients  between  master  and  slave  events  and  multiple  pass  culling 
of  the  rms  residuals  from  the  inversion  results.  This  effort  clearly  reduces  the  efficiency  with 
which  these  techniques  can  be  applied  to  large  clusters  and  presents  some  new  challenges 
in  developing  practical  algorithms.  However,  in  the  case  of  the  Larderello  1993  cluster,  we 
are  encouraged  that  the  final  result  so  clearly  agreed  with  know  geologic  structures.  This 
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suggests  that  refinements  of  the  correlation  techniques,  in  ways  that  reduce  the  effects  of 
scattering,  may  further  improve  the  results.  One  possible  direction  may  be  to  introduce 
more  thorough  modeling  of  the  scattered  wavefield  and  removal  of  the  non-essential  phases 
from  the  waveform.  This  may  improve  the  correlation  coefficients  by  increasing  the  usable 
window  length. 
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Figure  0-2:  The  Larderello  seismic  network  (triangles)  and  the  location  of  the  cluster  as 
determined  by  ENEL.  The  origin  is  located  at  Latitude  43N2’3T\  Longitude  10E29"5R". 
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Figure  0-4:  Five  synchronized  records  from  the  station  POMA  showing  P  and  S  arrivals. 
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Figure  0-5:  Cross-covariances  of  waveforms  from  five  events  with  the  waveform  from  event 
147  (type  II).  Note  that  the  record  for  147  is  therefore  the  auto-covariance  function. 
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Figure  0-6:  Cross-covariances  of  waveforms  from  five  events  with  the  waveform  from  event 
150  (type  I).  Note  that  the  record  for  150  is  therefore  the  auto-covariance  function. 
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Figure  0-7:  Amplitude  spectra  for  event  117  at  station  CORN. 
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Figure  0-8:  Master  events  of  Type  l  for  station  FOMA. 
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Figure  0-10:  P  wave  windows  for  a  master  and  slave  event  pair  (unfiltered). 
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Figure  0-12:  Auto-  and  cross-covariances  of  P  wave  windows  for  a  master  and  slave  pair 
(filtered  and  window-tapered). 
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Figure  0-13:  S  wave  windows  for  a  master  and  slave  event  pair  (unfiltered). 
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Figure  0-15:  Distribution  of  cross-correlation  offsets  for  P  wave  windows  (86  events  at 
station  POMA). 
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Figure  0-16:  Distribution  of  cross-correlation  offsets  for  S  wave  windows  (86  events  at  station 
POMA). 
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Figure  0-17:  Map  view  of  86  relocated  events  from  the  March  Larderello  microearthquake 
cluster. 
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Figure  0-18:  Local  map  view  of*  the  SO  relocated  events. 
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Figure  0-19:  Cross-sectional  view  of  the  86  relocated  events,  view  looking  west  in  depth 
section. 
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Figure  0-20:  Cross-sectional  view  of  the  86  relocated  events,  view  looking  north  in  depth 
section. 
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Figure  0-21:  Distribution  of  traveltime  residuals  from  the  relative  event  location  inversion 
of  86  events. 
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Figure  0-22:  Map  view  of  events  in  sub-set  1,  events  116  through  140. 
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Figure  0-23:  Map  view  of  events  in  sub-set  2,  events  141  through  156. 
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Figure  0-25:  Map  view  of  events  in  sub-set  4,  events  188  through  218. 
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Figure  0-26:  Estimated  error  ellipses  for  selected  events  from  subset  157-1 
values. 
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Figure  0-27:  Estimated  error  ellipses  for  selected  events  from  subset  188-218  showing  typical 
values. 


Of 


November  20,  1.997  19:01 


Figure  0-28:  Estimated  error  bars  for  selected  events  from  subset  157-18/  showing  typical 
values. 
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Figure  0-29:  Estimated  error  bars  for  selected  events  from  subset  188-218  showing  typical 
values. 
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Figure  0-30:  Focal  mechanisms  for  four  events.  These  examples  are  representative  of  t  he 
range  of  solutions  found  for  seventeen  events  from  the  March  1993  cluster.  Mechanisms  for 
events  117,  181,  and  187  are  average  solutions,  where  as  the  mechanism  for  event  200  are 
three  solutions  which  fit  the  observed  polarities  for  that  event. 
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Figure  0-31:  Epicentral  locations  of  the  event  subset  157-187  in  relation  to  depth  contours 
fo  the  seismic  “I<”  horizon  and  the  sub- vertical  discontinuities  (countour  map  source  ENEL. 
1988). 


Distance  Nortn  of  Origin  (km) 


28 


®  Pomarance 


32  33  34  35  36 

Distance  East  of  Origin  (km) 


Figure  0-32:  Epicentral  locations  of  the  event  subset  188-218  in  relation  to  depth  contour; 
fo  the  seismic  “K”  horizon  and  the  sub-vertical  discontinuities  (countour  map  source  KXFL 
1988). 
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